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Abstract 

The Lagrangian- Averaged Navier-Stokes alpha (LANS-a) model is a turbulence pa- 
rameterization that has been shown to capture some of the most important features 
of high resolution ocean modeling at lower resolution. Simulations using LANS-a in 
the POP primitive-equation ocean model resemble doubled- resolution simulations 
of standard POP in statistics like kinetic energy, eddy kinetic energy, and potential 
temperature fields. The computational cost of adding LANS-a is only 27% for our 
most efficient implementation, as compared to a factor of 8 — 10 for a doubling of 
resolution. 

The LANS-a model improves turbulence statistics with an additional nonlinear 
term and a smoothed advecting velocity. In this work we investigate different kinds of 
smoothing techniques and their effect on the LANS-a model's results and efficiency. 
We show that we can substitute convolution filters for full Hemlholtz inversions and 
produce similar results at a significantly lower expense. 

When constructing filters for LANS-a in a primitive-equation ocean model the 
filter weights must be chosen carefully, otherwise a pressure- velocity instability will 
be excited. We show analytically that certain ranges of filter weights are unstable, 
and confirm this with numerical experiments. Our stability criterion also guarantees 
that the kinetic energy is well defined, and that the filtered velocity is smoother that 
the original velocity. 



1 Introduction 



The LANS-a model is a turbulence model that can reproduce some of the statistics 
and structures that are seen in standard simulations (without LANS-a) at higher 
resolution. The LANS-a model accomplishes this in part by using a smoothed velocity 
field for the advecting velocity in the momentum and tracer equations. The question 
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addressed in this paper is, "What is the best way to smooth the velocity field in LANS- 
a?" This is not simply a matter of choosing the smoothing method that produces the 
best result. As a turbulence model, LANS-a must be extremely efficient as well; if the 
addition of LANS-a were expensive, it would be more sensible to drop the turbulence 
model altogether and simply run at higher resolution. Thus the choices we make in 
designing the LANS-a algorithm and smoothing filters are governed by the trade- 
offs of higher-resolution effects versus the cost of adding LANS-a to the model. We 
win this game if we can get the effects of, say, a doubling of resolution while only 
increasing the computation time by a small fraction of the time required to double 
the resolution with the standard model. 

We have designed a LANS-a algorithm and smoothing filter for primitive-equation 
ocean models that meets this criteria of a successful turbulence model. The algorithm 
is presented by Hecht, Holm, Petersen, and Wingate in a closely related work [1]. The 
focus of the present paper is the method used to smooth the advecting velocity. In 
the typical LANS-a model derivation, the smooth velocity, u, is computed from the 
rough velocity, v, with an inverse Helmholtz operator, 



where the namesake a parameter is a length scale that determines the amount of 
smoothing [2] . The Helmholtz inversion is the basis of comparison for other smoothing 
methods in this work, as this is the relationship between u and v in the original LANS- 
a equations. With it we are able to run successful simulations of LANS-a; however, 
it requires an iterative solver and is therefore relatively expensive. Fortunately, local 
averaging methods work well in LANS-a and are much cheaper. In these methods, 
the smooth velocity is computed as a weighted average of nearby neighbors, and can 
be written as a convolution with the filter function g as 



The LANS-a equation using a filter (sometimes called the Kelvin-filtered Navier- 
Stokes equation) retains Kelvin's circulation theorem [3], which is a property of funda- 
mental importance to LANS-a. If the filter is properly designed, other characteristics 
of LANS-a with the Helmholtz inversion are guaranteed: u is smoother that v, the 
energy is well defined, and LANS-a conserves energy in the absence of dissipation. In 
this paper we develop design criteria for filters that satisfy these requirements. 

There is a precedent of using simple filters in LANS-a in large eddy simulation (LES) 
models. Geurts and Holm [4] used a simple three-point top-hat filter as a smoothing 
operator in LANS-a and Leray LES modeling of three-dimensional isotropic turbu- 
lent mixing. They found that the LANS-a and Leray models are considerably more 
accurate than dynamic eddy- viscosity models, and at a lower computational cost. The 
LANS-a and Leray models performed particularly well at capturing the flow features 
characteristic of the smaller resolved scales. 

This paper is organized as follows: §2 briefly reviews the LANS-a equations for 
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primitive-equations and the POP-a algorithm. The pressure- velocity instability, which 
develops with the incorrect choice of filter weights, is discussed in §3. Using an eigen- 
value analysis of the discretized filter in §4, we develop design criteria for the filter 
weights to ensure the proper smoothing and kinetic energy properties. The results 
section, §5, shows that POP-a can produce results qualitatively similar to those of 
higher resolution simulations of standard POP using either the Helmholtz inversion 
or an filter to smooth the velocity. In the final section, §6, we compare the merits of 
each. The family of filters presented here produce stable simulations and are 20 — 50% 
faster than Helmholtz inversions. Filters produce better temperature statistics in the 
baroclinic instability model problem, while the Helmholtz inversion produces higher 
kinetic and eddy kinetic energy. 



2 The LANS-o model 



This section gives a brief overview of the primitive-equation version of LANS-a that 
is used in the POP-a code. These equations were originally derived by Holm et al. 
[5], and are described at length, along with the POP-a algorithm, by Hecht et al. [1]. 
The model equations are 
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where (vi, v 2 , v 3 ) = (v, v 3 ) and (u\, u 2 , u 3 ) = (u, u 3 ) are the rough and smooth veloc- 
ities, V is the horizontal gradient, <p is a tracer (temperature, salinity, or a passive 
tracer), p is the full density, po is the background density, n is a modified pressure, 
p is the physical pressure, T and T> are momentum and tracer diffusion terms [6], f 
is the Coriolis parameter, g is gravitational acceleration, and a is the alpha model's 
smoothing length scale. These equations are hydrostatic (4), incompressible (6), and 
Boussinesq (p appears in eqn 3 but not in 4). 

The principle feature of the LANS-a model is that the momentum equation (3) is an 
advection-diffusion equation for a Lagrangian-averaged velocity v, while the advecting 
velocity is an Eulerian-averaged velocity u. These names arise in the derivation of the 
LANS-a model [2], where velocities are averaged along a particle track (Lagrangian) 
or at a particular location (Eulerian). The Helmholtz relation (7) indicates that u 
is smoother that v. In this paper we prefer the names rough velocity for v, which is 
computed prognostically from (3), and smooth velocity for u, which is computed diag- 
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nostically from (7). The LANS-a model has an additional nonlinear term Y^=i v j^ u j 
in the momentum equation (3), which does not appear in the Leray model [2], and 
which allows LANS-a to satisfy Kelvin's circulation theorem while conserving energy 
and a form of potential vorticity in the absence of dissipation [2]. 

In this paper, we investigate how to most effectively and efficiently smooth the rough 
velocity v to obtain the smooth velocity u. This could be with a Helmholtz inversion, 
as in (7), or with a convolution using various filter functions, as in (2). LANS-a using 
a filter rather than a Helmholtz inversion still satisfies Kelvin's circulation theorem 
[3]. If the proper filters are chosen, one also conserves energy and a form of potential 
vorticity in the absence of dissipation. 

The most commonly used form of the LANS-a equations use a three-dimensional (3D) 
smoothing. In the primitive equations, which assume a thin layer approximation, this 
reduces to a 2D smoothing, and so our filters are horizontal-only as well. 3D smoothing 
may be tested at a future time, but we expect that the computational cost would be 
too high to justify the benefit. The Helmholtz inversion is horizontally isotropic. Thus 
one of the design criteria for the smoothing filter is that it is horizontally isotropic as 
well. 

The POP primitive-equation ocean-climate model [7,8], developed and maintained at 
Los Alamos National Laboratory, was augmented to include the LANS-a formulation 
[1]. POP is a split implicit /explicit code; that is, it takes an implicit step for the 
barotropic (vertically integrated 2D) velocity field, and takes an explicit step for the 
baroclinic (remaining 3D) velocity field. This is done so that the timestep is not 
limited by the fast free-surface gravity waves in the barotropic velocity field. 

In the LANS-a version of POP, named POP-a, a smooth velocity field must be com- 
puted for both the baroclinic and barotropic velocity fields. In the explicit baroclinic 
component of the code, the smooth baroclinic velocity, u, is computed by simply 
smoothing v. The implicit barotropic component is more complicated, and in the 
end a reduced version of the full barotropic algorithm, as derived from the LANS-a 
equations, was used [1]. The reduced algorithm produced nearly identical results to 
the full algorithm, but was three to four times faster. 

A few pieces of the POP-a barotropic algorithm are needed to understand the pressure- 
velocity instability presented in the following section. The 2D barotropic velocities, 
denoted by capital letters, are found by integrating the full velocity from the bottom 
to the free surface height, rj: 



= — ^— r u dz, 

H + 77 J-H 



where H(x,y) is the ocean depth when the surface is at rest. By integrating the 
continuity equation (6), which states that the smooth velocity is divergence- free, we 
obtain a prognostic equation for the free surface height, 
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[ V (V ■u + d z u 3 )dz = 0, (10) 

J-H 

^L + V .(H + V )V = 0. (11) 

The free surface height now replaces the pressure in the barotropic form of the mo- 
mentum equation, 

™ + fxV = -gVri + G, (12) 

where G contains the vertically integrated forcing terms. The barotropic momentum 
equation is the same as the shallow water equation. 



3 Pressure-velocity instability 

In this section, we present an analytic stability analysis of a pressure-velocity insta- 
bility that constrains the filter weights in the barotropic filter. The pressure-velocity 
instability is a phenomenon particular to models where a smoothed velocity field is 
used in the equation for free surface height r\ (11). Consider a velocity field in one- 
dimensional (ID) discretized shallow water equations, where the pressure and velocity 
gridpoints are offset. If the velocity field were at the Nyquist frequency (alternating 
between positive and negative at each grid cell), then using the standard shallow 
water equations, the free surface height would raise (lower) where the velocity field 
converges (diverges). In either case, a pressure gradient force would ensue that would 
work against the tendency of the surface to undergo further deformation, providing 
a stabilizing effect. 

In the LANS-a model the equation for the free surface height involves the smooth 
velocity because this equation is derived from the continuity equation (6), which 
involves the smooth velocity. When the outer stencil weights are sufficiently large in 
LANS-a, a pressure-velocity instability ensues in the highest frequencies. This can be 
shown with the simplest ID velocity field and a smoothing stencil of width 

three, 



bvi-x + avi + bv i+1 

Ui = aTYb • (13) 

In most cases, Ui and V{ will have the same sign at each gridpoint, as one would expect 
for a smoothing operator. However, for high frequencies and sufficiently large outer 
stencil weight b, u can have the opposite sign of v, as shown in Fig. 1 for the Nyquist 
frequency and b = 0.7. 

The tendency equation for v depends on 77, but 77 in turn is not directly dependent 
on the convergence/divergence of v but on that of u. In the case of a pure harmonic 
mode, restricting the local filter to produce a smoothed velocity u of the same sign 
everywhere as the rough velocity v is sufficient to also ensure that the divergence 
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of smooth and rough velocities share the same sign, preserving the essential physical 
feedback of the underlying equations even after inclusion of the LANS-a model. Cases 
in which this physical feedback is respected or not, depending on the filter width, are 
illustrated in Fig. 2. 

This understanding of the pressure-velocity instability allows us to construct an an- 
alytical stability criterion as follows: The filter weights must be chosen so that the 
filtered velocity u retains the same sign as v for all harmonic modes resolved on the 
grid. 

For example, for the stencil of width three (13) where a is normalized to 1, consider 
the Nyquist frequency 



v i+j = cos(nj/k) (14) 

with k — 1. Since Vi = 1, the condition that u and v are like-signed is Ui > 0. 
Using (13), the constraint on the stencil is b < 0.5. For a slightly higher wavenumber, 
k = 1.5, the constraint is weaker, b < 1. For still higher wavenumbers, u and v are 
like-signed for any choice of b. Thus b must be chosen based on the tightest constraint, 
b < 0.5. 

Numerical simulations of POP-a with a stencil of width three almost exactly cor- 
respond to the predictions of this stability analysis. The model runs stably when 
b < 0.48. At larger values of b an instability quickly grows in v and rj at the Nyquist 
frequency, so that v is four times larger than u within 100 time steps (Fig 3). By 
150 timesteps the conjugate gradient routine in the barotropic solver fails to con- 
verge within 1000 iterations. Numerical simulations use a 2D stencil that is a squared 
version of the ID stencil, as discussed below. 

Because of the limitation imposed by the pressure- velocity instability, further smooth- 
ing must be accomplished by increasing the stencil size. The criterion that u and 
v must be like-signed for all frequencies applies, but larger stencils have more un- 
knowns, one for each stencil weight. Figure 4 shows these constraints in chart form 
for ID filters of width 3, 5, 7, and 9 and for wavenumbers k = 1 ... 4. In general, 
these constraints require that the weights are each less than 0.5. Beyond that, each 
consecutive weight moving away from the center must be slightly less than the one 
before it. This latter requirement can be seen in the constraints for the stencil of 
width nine for wavenumber k = 3 and k = 4, 



c + 2d + e<l + b, (15) 

V2d + 2e < 1 + V2b. (16) 

Neither of these are satisfied if all the weights are 0.5, but are satisfied when the 
weights decrease as 



6 = 0.45, c = 0.4, d = 0.35, e = 0.3. (17) 

Again, the results of this stability analysis were verified using the POP-a model. 
Simulations using the stencil weights in (17) were always stable (Fig 5). If any of the 
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weights are increased by 0.02, a fast-growing, high frequency appears in v and rj (Fig 



In the POP-a model the smoothing operator is 2D. The 2D stencil was chosen to be a 
squared version of the ID stencil (Fig. 6). Such a stencil is isotropic, and thus follows 
the same stability analysis as the ID case presented above. POP-a is very sensitive 
to these stability constraints; for example, if the corner weights of the stencil in Fig. 
6b are increased by 10 — 20%, the simulation is unstable. 

Other experiments were conducted where the 2D stencil was a diamond or circle 
instead of a square. The motivation here is that in the standard square 2D stencil 
(Fig. 6) the corner weights are small, and don't contribute much to the average. Thus 
computation time could be saved by removing the corners from the stencil, and the 
degree of smoothing should be nearly the same. However, in practice any stencil other 
than a full square stencil resulted in unstable simulations that did not converge in the 
barotropic solver. This was even true when only a single gridpoint was left off of the 
stencil on each corner, for example when the 0.09 weight is removed from a nine-wide 
stencil in Fig. 6b. Thus all of the simulations presented in this work use a complete 
square stencil. 

Near boundaries, the stencil is reduced in size so that all grid-cells used in the aver- 
aging calculation are water-cells, and none are land-cells. For example, at a grid-cell 
directly beside or diagonal to a land cell no averaging takes place (i.e. Uij = Vij). 
If the nearest land-cell is two grid-cells away, a filter of width three is used. The 
disadvantage of this method is that the smoothing goes to zero as one approaches a 
land-cell. This stands in contrast to the LANS-a equations, where the smoothing scale 
a is constant throughout the domain. Other filtering schemes near the boundary were 
considered: One could use a constant-size stencil and assign zero velocity to land-cells, 
but this strongly damps the smooth velocity near the boundary. Another possibility 
is to use large but nonsymmetric stencils near the boundary, but this proved to be 
unstable. 

4 Eigenvalue analysis of the filter 

Most theoretical results for the LANS-a model depend on the two velocities being 
related through the Helmholtz operator. We can replace the Helmholtz operator with 
filters if two critical properties are preserved. These are that (1) the filtered velocity 
u must be smoother than the unfiltered velocity v, and (2) that the energy, 



must be well defined (e.g., non-negative). When both of these criterion are satisfied, 
the filter may be used in place of the Helmholtz inversion operator in the LANS- 
at model. In this section we will review why the Helmholtz inversion satisfies these 
criterion, and then proceed to show that the filter described in the previous section 
satisfies them as well. 



3). 




(18) 
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The Helmholtz inversion is easily shown to be a smoothing operator using Fourier 
analysis, 



u(k,l) 



v(k,l) 



a 2 k 2 ' 



(19) 



where u and v are the Fourier coefficients, (k, I) are horizontal wavenumbers and 
k\ = k 2 + I 2 . For scales much larger than a the Fourier coefficients are unaffected, 
because 1/kh » a so a 2 /^ << 1 and u ~ v. For scales at or smaller than a the 
Fourier coefficients are strongly damped because a 2 k 2 > 1 so that u < v. The energy 
for the Helmholtz inversion is also well defined, because (18) can be rewritten, using 
integration by parts, as 



E= l - 
2 



2 . 2 |V7 1 2 t3 

+ a Vu ax, 



(20) 



which is nonnegative. 



We next show that these properties also hold for our filters by using the eigenvalue 
problem applied to the discrete filter-operator matrix A. To illustrate this process, 
first consider the simple three-wide filter in ID shown in (13). This filter can be 
written as a tridiagonal matrix, 



A 



1 



a + 2b 



a b ■•• 

b a b ••■ 

o '•• '•• ■■. 

■■• b a 



(21) 



so that a vector can be filtered using u = Av. The eigenvalues of A vary as a function 
of b. If b = 0, A is the identity matrix and the eigenvalues are all one. As b increases, 
the range of eigenvalues increases but always remains less than one (Fig. 7). The 
smallest eigenvalues are associated with the highest wavenumber modes (Fig. 8), 
indicating that the most oscillatory modes are damped the strongest, and that the 
matrix A is indeed a smoothing operator. 

We now evaluate the eigenvalues of the full filter (Fig. 6) in 2D. On an n x n grid 
the filter operator is an n 2 matrix (Fig. 9). Here we consider Dirichlet boundary 
conditions of zero, but the results are similar for periodic boundaries. As in the 
ID case, all eigenvalues are less than one. We have chosen weights so that most 
eigenvalues are near zero but positive (Fig. 10). The most oscillatory modes have 
very small eigenvalues, so that these are strongly damped (Fig. 11). This ensures that 
the filter is a smoothing operator. 

The last task is to ensure that the energy defined in (18) is well defined. The matrix 
A is positive definite if v T Av > for all nonzero vectors v. Since the filtered velocity 
u = Av, this is equivalent to v • u > 0, which means that the global energy (18) 
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is nonnegative and therefore well defined. A matrix is positive definite if all of its 
eigenvalues are positive. Thus a criterion for our filter is that its minimum eigenvalue 
is positive. This leads to a process for choosing each filter weight: plot the minimum 
eigenvalue as a function of the filter weight, and choose the largest weight that guar- 
antees stability. There is no reason to choose a smaller value, because then the filter 
smooths less. Figure 12 shows that this process produces the same restrictions on 
filter weights as the pressure- velocity instability in the previous section. In fact, the 
two methods of analysis agree so well that the pressure-velocity instability appears 
to be a physical manifestation of the poorly defined energy. Numerical experiments 
where the filter weights were varied confirm that the simulations are unstable when 
the energy is poorly defined. 



5 Results 

The model problem is a zonally periodic channel with a deep-sea ridge, eastward 
wind forcing, and surface thermal forcing that restores SST to a smooth profile from 
12°C in the north to 2°C in the south, as described in [1]. This configuration is 
an idealization of the Antarctic Circulpolar Current, and was designed to induce 
baroclinic instability, where isopycnals are tilted from the horizontal meridionally 
by the surface thermal and wind forcing. If present, eddies in the flow advect heat 
meridionally, with the net effect that the isopycnals are less tilted; in this way, the 
eddies convert the potential energy of the tilted isopycnals to the kinetic energy of 
the eddies themselves. Simulations of standard POP at three resolutions show this 
effect. In the lowest resolution simulation (0.8, see Table 1) the Rossby Radius of 
deformation is not resolved, and no eddies are present. As the resolution doubles 
(0.4) and doubles again (0.2) more eddies exist, and the isopycnals are less tilted 



The flatter isotherms at higher resolution create a sharper thermocline and colder 
temperatures below the thermocline (Fig. 14). (In this model configuration, salinity 
is constant, so the isotherms shown in the figures coincide with isopycnals.) Higher 
resolution simulations of standard POP also show progressively higher kinetic energy 
(Fig. 15) and higher eddy kinetic energy (Fig. 16) as more eddies are resolved. 

POP-a simulations resemble higher resolution simulations of standard POP in all of 
these measures. The strength of the eddy activity is controlled by how much smooth- 
ing is done on the rough velocity v: when a Helmholtz inversion is used, a larger 
a smooths more; when a filter is used, wider filter stencils smooth more. In either 
case, Figs. 13-16 show that POP-a can produce results that look progressively like a 
doubling of resolution of standard POP as the smoothing is increased. 

Because POP-a has two velocities, some details are necessary to explain how these 
diagnostic quantities were computed. The kinetic energy is 



(Fig. 13). 




(22) 
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for POP and POP-a, where the subscripts indicate horizontal components. The eddy 
kinetic energy is 



EKE = («) 2 + K) 2 ) /2, EKE a = (u' lV 'i + u' 2 v' 2 ) 1% (23) 

where each variable is the sum of a mean and perturbation, e.g. U\ — ul + u' x , and 
the overbar indicates a five-year time-averaged mean. 

There would seem to be three possible ways one could compute the kinetic energy 
for POP-a: using the products u 2 , uv, or v 2 . The product of the smooth and rough 
velocities, uv, is used because this the conserved kinetic energy that emerges from 
the derivation of the LANS-a equation [5,2]. For comparison, kinetic and eddy kinetic 
energies were also computed using u 2 and v 2 . The product of the smooth velocities is 
smaller than uv, and the product of the rough velocities is larger (typically by 10 to 
20%). However, the general trends shown in Figs. 15 and 16 still hold if u 2 or v 2 are 
used instead of uv. 

An interesting difference between the filter and Helmholtz inversion in POP-a is that 
simulations using filters are better at the temperature statistics, while simulations 
using the Helmholtz inversion have higher kinetic and eddy kinetic energy. This can 
be seen by comparing 0.8H2 with 0.8F9 or 0.4H1.5 with 0.4F9 if figures 13 — 16. 
Even though results differ a bit with each measure, one may loosely assign each 
filter an effective alpha value. For example, in most plots the filter F5 corresponds to 
a = 1.5Ax and F9 corresponds to a = 2Ax. 

The computation time required for POP-a is slightly longer than standard POP due 
to the smoothing operations and additional nonlinear term, but is much less than a 
doubling of resolution in standard POP (Fig. 17). Computation increases as the filter 
width increases from three to nine; smoothing using the Helmholtz inversion is even 
slower than the nine-wide filter. 

The standard POP algorithm may use an explicit or implicit discretization for the 
Coriolis term. Typically the implicit discretization is chosen because it allows longer 
timesteps to be taken, resulting in faster simulations. As described in [1], the POP- 
a algorithm must use the explicit implementation of the Coriolis term. Fortunately, 
POP-a can use longer timesteps than the standard POP algorithm with explicit 
Coriolis [1,9]. This explains why most of the POP-a simulations in Fig. 17 are faster 
than POP with an explicit Coriolis term, and just a bit slower than POP with an 
implicit Coriolis term. 

6 Conclusions 

This paper presents an assessment of two methods of smoothing the velocity field 
in the LANS-a turbulence model: the Helmholtz inversion and filters. The LANS-a 
model equations specify the Helmholtz inversion, but this method requires an iterative 
conjugate gradient routine for each smoothing. The filter, which is simply a weighted 
average of nearby neighbors, results in simulations that are 20 to 50% faster than 
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those with the Helmholtz inversion. This disparity is expected to increase with larger 
domains and more processors, because more iterations are required for convergence 
in the conjugate gradient method as the problem-size increases ([10], p. 530). The 
savings in computation time provides a strong motivation to choose a local filter over 
the global Helmholtz inversion. 

Both the Helmholtz inversion and the filters produce results expected of the LANS-a 
turbulence model: statistics such as temperature profiles, kinetic and eddy kinetic en- 
ergy resemble higher-resolution simulations of non-alpha simulations. The parameter 
a is a length-scale that controls the amount of smoothing, and thus the strength of 
the turbulence model. Analogously, as the filter width varies from three to nine, the 
smoothing increases and the effects of LANS-a are stronger. For a particular value of 
a or filter width, the methods perform a bit differently in different metrics: the filter 
has better temperature profiles (as judged by higher-resolution simulations), while 
the Helmholtz inversion has higher kinetic and eddy kinetic energy. 

The filter is more robust than the Helmholtz inversion. When a was too high (i.e the 
smoothing too strong) in the Helmholtz inversion method, the kinetic energy grew 
without bound and the simulation was unstable. This is typical of LANS-a models. 
Unstable simulations occurred when a = 2.5Ax and 2Ax at resolutions of 0.8 and 0.4, 
respectively. For the filter, all stencils tested (up to 9x9) were stable, and stronger 
sub-grid model effects could be obtained with the filter than with the Helmholtz 
inversion. 

The stencil weights must be chosen carefully, otherwise POP-a will be unstable. If 
the outer weights are too large, the smooth velocity will be of the opposite sign as the 
rough velocity for the highest wave-number modes. When this occurs the pressure 
gradient, which is calculated from the smooth velocity, is of the wrong sign. The 
physically correct relationship between pressure gradient force and rough velocity 
is lost, and the mode grows. The conditions that bring about this pressure- velocity 
instability were identified exactly, and then verified by numerical simulations. We 
found that the first neighboring weight must be less than one-half the central weight, 
and each consecutive neighbor must be less than the one before. In 2D, this pattern 
is squared to make an isotropic, stable stencil. The full square was required; leaving 
off the corners to make diamond-shaped stencils resulted in unstable simulations. 

It is important that the filters chosen are actually smoothing operators, and that 
the energy is well defined. Both of these properties were verified using an eigenvalue 
analysis of the filter, discretized as a matrix operator. In fact, the criterion on filter 
weights that guaranteed well-defined energy turned out to be the same as the weight 
limits to prevent the pressure- velocity instability. Thus it appears that the instability 
is a physical manifestation of poorly-defined energy in the equation set. 

The goal of any turbulence model is to capture the effects of higher resolution simula- 
tions without paying the computational price of running at that higher resolution. We 
have shown that our implementation of the LANS-a model in a primitive-equation 
ocean-climate model accomplishes this goal. Specifically, statistics in POP-a simula- 
tions such as kinetic energy, eddy kinetic energy, and temperature profiles are similar 
to standard POP simulations with double the resolution. But the LANS-a simula- 
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tions only required 27% more computing time than standard POP, as compared to a 
factor of nine increase to double the resolution of standard POP (for the 0.4 case). 

We are currently performing simulations to compare LANS-a with other turbu- 
lence models, including constant coefficient hyperviscosity, Gent-Mc Williams isopyc- 
nal tracer mixing [11], the Leray model [2], and the Simplified Bardina model [12]. 
These results, to be published in the near future, will consider the merits of LANS-a 
relative to competing models in the context of primitive-equation ocean models. We 
also intend to test the model in an ocean-basin or global simulation. 

The comparison of the Helmholtz inversion and the filter in this study has shown 
that the filter is both cheaper and more robust than the Helmholtz inversion, while 
producing similar results. Thus our future work with the POP-a model will use local 
filters in place of the global Helmholtz for the majority of our simulations. 
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Fig. 1. The effects of filtering the rough velocity v (open circles) to compute the smooth 
velocity u (closed circles) for a ID filter of width three, where the weight of the neighboring 
cell, 6, is varied between 0.3 and 0.7, as shown in the left column. For the Nyquist frequency 
(k = 1, middle column) u is of opposite sign to v if b > 0.5; this results in an instability, 
since the equation for free surface height i] involves u. For a slightly higher wavenumber 
[k = 1.5, right column) u has the same sign as v for b < 1. 
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Fig. 2. Schematic of variables to illustrate the pressure-velocity feedback instability. In the 
stable case (a) the smooth and rough velocities, u and v, are of the same sign, and the 
pressure gradient force — Vr? counters the original velocity perturbation, as expected. In the 
unstable case (b), u and v have opposite signs; the free surface height, which is computed 
from u, increases (decreases) where u converges (diverges); now the pressure gradient force 
reinforces the rough velocity v , increasing this perturbation. 
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Fig. 7. Eigenvalues of the ID, width-three filter operator using 16 gridpoints. When 6 = 
the operator is simply the identity, so all eigenvalues are one. As b increases the eigenvalues 
spread out below one, indicating that some eigenvectors contract. Eigenvalues are never 
larger than one, so no eigenvectors increase. When b > 0.5 there are negative eigenvalues. 




Fig. 8. Eigenvectors of the ID, width-three filter operator using 32 gridpoints. High 
wavenumber modes (a,b) have eigenvalues near zero, while low wavenumber modes (c,d) 
have eigenvalues near one. This shows that the filter is a smoothing operator, because the 
highest wavenumbers are strongly damped. 
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Fig. 9. Matrix representation of the 2D, 7x7 filter using 10x10 gridpoints. The full matrix 
is 100 elements square, while each block is 10 square. This matrix is not yet normalized by 
the factor l/(a + 26 + 2c + 2d + 2e) 2 , so weights range between and 1. 
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Fig. 10. Distribution of eigenvalues of the 2D, 9x9 filter operator using 16x16 gridpoints. 
Filter weights are as shown in Fig. 5. The majority of eigenvalues are less than 0.1, indicating 
that most eigenmodes are stongly damped. 
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Fig. 11. Eigenmodes of the 2D, 9x9 filter operator using 16x16 gridpoints. As in the ID 
case, highly oscillatory modes (a) have eigenvalues near zero and so are strongly damped, 
while less oscillatory modes (d) have eigenvalues near one and so are only weakly damped. 
This shows that the filter we have implemented is a smoothing operator. 
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Fig. 12. Minimum eigenvalue of the 2D, 9x9 filter operator using 16x16 gridpoints, as a 
function of each filter weight. If all eigenvalues are strictly positive then the operator is 
positive definite, and the energy is well defined. The value of each filter weight is chosen 
(starred points) to satisfy this criterion. 




Fig. 13. Depth of 6°C isotherm of potential temperature using the Helmholtz inversion (a) 
and the filter (b) to smooth the velocity. Higher resolution simulations of standard POP 
(dotted line) have isotherms which are less tilted than lower resolution simulations (solid 
lines). As a or filter size increases, the POP-a simulations show flatter isotherms, and 
approach a doubling of resolution using standard POP. 
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(a) Helmholtz inversion (b) filters 

Fig. 14. Horizontally averaged potential temperature versus depth using the Helmholtz 
inversion (a) and the filter (b) to smooth the velocity. Higher resolution simulations of 
standard POP (dotted line) are cooler and have a sharper thermocline than lower resolution 
simulations (solid lines). Again, POP-a simulations approach a doubling of resolution using 
standard POP. 
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Fig. 15. Kinetic energy for all simulations. As the resolution increases with standard POP, 
the kinetic energy increases. Kinetic energy also increases using POP-a at fixed resolution, 
when a is increased using the Helmholtz inversion (a) or the stencil width of the filter is 
increased (b). Each value was calculated by averaging over the entire domain, and in time 
after 500 years. 
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Fig. 16. Same as Fig. 15 but for eddy kinetic energy. Note that at the lowest resolution 
the eddy kinetic energy of standard POP is indistinguishable from zero on this scale. 
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Fig. 17. Timing data from various implementations of POP and POP-a (note log scale). 
The POP-a filter algorithms (F3 — F9) are much cheaper than a doubling of resulution using 
standard POP. Using a Helmholtz inversion to smooth is more expensive than the filters. 
All POP-a algorithms use explicit Coriolis discretization. 
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Table 1 

Model parameters for experiments discussed in this paper, fw is the filter width; grid is 
the number of gridpoints in (x,y,z); Ion is the longitudinal grid-cell width; and lat is the 
latitudinal grid-cell width. The names are a concatenation of the meridional resolution, the 
type of smoothing, and the value of a or filter width. 
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Table 2 

Minimimum steps/day and the resulting clock time for various algorithms and resolu- 
tions. The second column states whether the barotropic Coriolis term is implicit or exlpicit. 
Clock time is in processor-hours per simulated decade. The fastest POP-a algorithm is the 
reduced algorithm with a filter. Even though the POP-a algorithms use explicit Coriolis 
discretization, the timestep is larger (fewer steps per day) than standard POP with explicit 
Coriolis. 
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